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ABSTRACT 

We perform a detailed comparison of the phase-space density traced by the particle 
distribution in Gadget simulations to the result obtained with a spherical Vlasov solver 
using the splitting algorithm. The systems considered are apodized Henon spheres with 
two values of the virial ratio, R ^ 0.1 and 0.5. After checking that spherical symmetry 
is well preserved by the A-body simulations, visual and quantitative comparisons are 
performed. In particular we introduce new statistics, correlators and entropic estima¬ 
tors, based on the likelihood of whether A-body simulations actually trace randomly 
the Vlasov phase-space density. When taking into account the limits of both the A- 
body and the Vlasov codes, namely collective effects due to the particle shot noise in 
the hrst case and diffusion and possible nonlinear instabilities due to hnite resolution 
of the phase-space grid in the second case, we hnd a spectacular agreement between 
both methods, even in regions of phase-space where nontrivial physical instabilities 
develop. However, in the colder case, R = 0.1, it was not possible to prove actual 
numerical convergence of the A-body results after a number of dynamical times, even 
with A = 10^ particles. 
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1 INTRODUCTION 


Stars in galaxies and dark matter in the Universe can be 
modeled in phase-space as self-gravitating collisionless fluids 
obeying the Vlasov-Poisson equations: 


%+U.Vrf-Vr(|>■V^.f = 0, 
at 


Arcf = AttGp = 47rC 


jtir, 


u,t) dit, 


( 1 ) 

( 2 ) 


where f(r,u,t) represents the phase-space density at posi¬ 
tion r and velocity it, cf is the gravitational potential, and 
G is the gravitational constant. 

In general, these equations do not have simple ana¬ 
lytical solutions. They are therefore often solved numeri¬ 
cally. The most widely used numerical scheme is the A- 
body approach and there exist many different implementa¬ 
tions, which mainly differ from each other in the way Pois¬ 
son equation is solved (see, e.g.,|Bertschinger|1998 iColombi 


|200T1 IDolag et al.||2008l [Dehnen fc Read||20il[ for reviews 

on the subject). The A-body method attempts to sample 
the phase-space density by an ensemble of Dirac functions 
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that represent particles interacting with each other through 
gravitational force. In order to avoid numerical artefacts due 
to the 1/r^ divergence of the force at small distances, the 
gravitational potential is usually replaced by an effective one 
so that the force is smoothed at scales smaller than a soft¬ 
ening parameter e. This procedure corresponds to assuming 
that the particles are clouds of size e interacting with each 
other. 


Approximating the phase-space density with macro¬ 
particles, however, has its own limitation. In particular, the 
close A-body encounter is one of the most notable sources of 
numerical artefacts, in addition to more subtle collective ef¬ 
fects induced by the discrete nature of the distribution of the 


particles (see, e.g. 

Aarseth, Lin, & Papaloizou||1988| 

|Splin- 

ter et al.|19981|Boi] 

y, Athanassoula, & Kroupa|2002||Binney 

2004| Joyce, Marcos, & Sylos Labini||2009|). Of course, the 


time integration scheme and the way to solve the Poisson 
equation numerically are well-known sources of errors, even 
though not particular to the A-body method. 

There are several previous studies that discussed the 
limitations of the A-body results, including underestimat¬ 
ing strong numerical artefacts, particularly in the cold case 


where the initial velocity dispersion is null (see, e.g., Melott 
et al.|1997|rMelott|2007| l, and long-term nonlinear resonant 
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modes induced by the discrete nature of the particles (see, 
e.g., Alard Colombi||2QQ5 Colombi Toum^|2014 ). We 
also note that it is not yet obvious that the fine inner struc¬ 
ture of dark matter halos is completely understood from 
physical and even numerical points of view, despite numer¬ 
ous intensive convergences studies of the Wbody approach 
(see, e.g., Moore et al.|199^ Jing Suto|20'00 2002 Power 


et al. 1981) and spherical systems (Fujiwara 1983). Nev- 


|et al.||2003 Springel et al.|2608 Stadel et al.||2609). 

It is therefore highly desired to develop alternative nu¬ 
merical methods to the traditional Wbody approach so that 
one can understand better its validity and fundamental lim¬ 
itations. 

In the cold case, relevant to the current paradigm of cold 
dark matter scenario, the phase-space distribution function 
is supported by a three-dimensional sheet evolving in six¬ 
dimensional phase-space, which can be partitioned in a con¬ 
tinuous way with an ensemble of tetrahedra as proposed in 
recent works (see, e.g., Shandarin, Habib, Heitmann|20T2 


Hahn, Abel, Kaehler|2013 ). Unfortunately, the increasing 

complexity of the structure of the system during evolution 
requires more and more sampling elements, and the com¬ 
putational cost becomes prohibitive after several dynamical 
time-scales. 

In this article, we consider the warm case, in which the 
system presents a non-negligible initial local velocity disper¬ 
sion component relative to gravitational potential energy. 
In this case, the phase-space distribution function has to be 
sampled on a 6-dimensional mesh, which makes again the 
computational cost very high. Therefore, we shall restrict 
to spherical systems, hence reducing the actual number of 
dimensions of the dynamical setup to three. 

There exist many methods to solve the Vlasov-Poisson 
equations in the warm case, mainly developed in plasma 
physics. One of the most famous solvers is the splitting 


algorithm of Cheng V Knorr (1976) and its numerous ex- 

tensions (see, e.g.|Shoucri V Gagne|1978| iSonnendriicker et 

al.||1999 |Filbet, Sonnendriicker, & Bertrand 

|2001| IBesse & 

Sonnendriicker |20031 Alard & Colombi||2005 

Umeda||2008| 

Besse et al. |2008[ |Crouseilles, Y 

ehrenberger, & Sonnen-] 

driicker|2010||Campos Pinto|2011| 

Rossmanith & Seal|201lf 

Giiglii, Christlieb, & Hitchon||2014 

but this list is far from 


complete). This algorithm, that we shall adopt below, ex¬ 
ploits directly the Liouville theorem: the phase-space density 
f{r,v,t) is conserved along motion. Then the equations of 
the dynamics during each time step are divided into “drift” 
and “kick” parts according to Hamiltonian dynamics and are 
solved backwards: 


rir,u) 

= f{r - 

uAtl2^ It, t), 

Drift, 

(3) 

r*{r,u) 

= Pir, 

U + Vr^At), 

Kick, 

(4) 

f{r,u,t + At) 

= PPr 

— itAt/2, It), 

Drift, 

(5) 


where Vt» 0 is computed from /*. In practice the phase-space 
distribution function is sampled on a mesh, and each step is 
performed by using tracer particles located at mesh sites and 
following the equations of motion split as above. Resampling 
of /*, /** and finally the phase-space distribution function 
at the next time step is performed by using an interpolation, 
e.g. based on the spline method. 

The splitting scheme was applied for the first time in 


astronomy in early I980’s, to one dimensional systems ( Fu- 
jiwarajpTM ), galactic disks ( Watanabe et al.||I98I Nishida 


ertheless, it has been almost forgotten since then except 
for a few contributions (e.g., Hozumi, Fujiwara, V Kan-Ya 


1996 Hozumi, Burkert, & Fujiwara 2000) that include a 


recent preliminary investigation of the algorithm in full 6- 


dimensional phase-space (Yoshikawa, Yoshida, V Umemura 


2013 ). 

As mentioned above, however, solving fully six¬ 
dimensional phase-space problems with sufficient accuracy 
is still very unrealistic now. In this article, therefore, we fo¬ 
cus on spherical systems, where phase-space is only three 
dimensional: the three coordinates of interest are the ra¬ 
dial position r, the radial velocity v and the angular mo¬ 
mentum j. Following earlier works performed in the frame¬ 
work of one dimensional gravity (see, e.g., Mineau, Feix, V 
Rouet 1990), we carry out a detailed comparison between 
an V-body code. Gadget ( Springel, Yoshida, V White|200I 


,Springel||2QQ5 ), and an improved version of the splitting al¬ 
gorithm implementation by Fujiwara (1983), VlaSolvelj 

Our goal is to check how well the particle distribu¬ 
tion in Gadget traces the phase-space density obtained from 
VlaSolve, and to see how the results depend on various pa¬ 
rameters of the simulations, in particular the number of par¬ 
ticles in the V-body simulations and the spatial resolution 
in the Vlasov code. We would however like to emphasize 
here that the purpose of this article is not to compare the 
performance of the two codes from the view-point of com¬ 
putational cost. 

While a fairly good physical insight is obtained through 
visual inspection of the resulting phase-space density plots, 
we also present a more quantitative comparison. To do so, 
we introduce correlators and entropic estimators based on 
a likelihood approach, ans ask whether the iV-body simula¬ 
tions can be considered as local Poisson realizations of the 
Vlasov code phase-space density. 

Because of our restrictive choice of the geometry of the 
system, it is important to simulate spherical configurations 
that are known to be stable against small anisotropic pertur¬ 
bations induced by the shot noise of the particles. Indeed, 
we shall use the public treecode Gadget without any spe¬ 
cific modification to enforce spherical dynamics. Although 
an alternative approach consisting in enforcing pure radial 
dynamics in Gadget (see, e.g., Huss, Jain, fc Steinmetz|I999 ) 
may facilitate comparisons with the Vlasov code, we do not 
adopt this approach in order to avoid any possible subtle 
biases in the analyses. 

In this respect, the Henon sphere ( Henon||I964 ) is par¬ 
ticularly suited for our purpose since it is known to preserve 
well its spherical nature during the course of dynamics even 
when being simulated with a iV-body technique and, in par¬ 
ticular, it is not prone to radial orbit instability (see, e.g.. 


Albada||I982[ [Hozumi, Fujiwara, fc Kan-Ya||I996| |Roy| 


& Perez 2004[ [Barnes, Lanzel, fc Williams||20Q9| ). In this 
configuration, the initial phase-space distribution function 
is isotropic and Gaussian distributed in velocity space and 
given by 


h{r,vj) = 


po 


(27rcrg)3/2 


exp 


I v‘ 
'2~ 


+f/p 


^ VlaSolve can be downloaded from the following web page: 
www.vlasix.org. 
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3 


r ^ Rn, (6) 

with (47r/3)po-RH = the total mass of the system. In the 
simulations discussed in this article, we work in units where 
G = 1, and the initial radius of the Henon sphere and its 
total mass are chosen to be 

M = 1, i^H = 2, (7) 

which fixes ay in equation once the virial ratio is given. 

We shall consider “warm” and “cold” settings, which cor¬ 
respond to the initial virial ratio R = \2T/W\ = 5RucrylM 
of ~ 0.5 and ~ 0.1, respectively, where T and W are the total 
kinetic and potential energy of the system. The two classes 
of initial conditions exhibit distinct features, in particular 
concerning the metastable state to which the system re¬ 
laxes through phase mixing. The warm system builds a core- 


p{r) ~ r (see, e.g. __ __ 

|1982| . In contrast, the cold system develops a more concen¬ 


Henon 

1964 

Gott 

1973 

van Albada 


trated smaller core (see, e.g., van Albada|1982| Sylos Labini 


2012), but never reaches a strictly stationary regime because 


a significant fraction of the mass acquires positive energy 


and escapes from the system (see, e.g., van Albada 1982 
Joyce, Marcos, Sylos Labini|200^ [Sylos Labini||2012 ). 

This article is organized as follows. In §[^we describe 

our Vlasov solver, VlaSolve. Section [^provides information 
about the V-body runs and the parameters used in Gadget. 
In § [^ we check that the iV-body simulations stay indeed 
spherical during evolution. Section presents a visual in¬ 
spection of the phase-space density, which is followed by a 
quantitative statistical analysis in § [^ Finally, §[^ summa¬ 
rizes and discusses our present results. 


2 THE VLASOV CODE: VLASOLVE 

Under spherical symmetry, the Vlasov equation reads 




dt 


dr 


dv 


( 8 ) 


where v is the radial component of the velocity, j is the 
angular momentum. My = M (< r) is the mass inside a 
sphere of radius r. 

Our code VlaSolve solves equation ^ numeric ally with 
the splitting algorithm, following closely Fujiwara (1983). 


Phase space is discretized into a rectangular mesh of 
size (iV-r , AT , 0 for -Rmin ^ “C ^ -Rmaxj "Cmax 'Cmax 5 

and 0 ^ j ^ e/max- More specifically, we use a logarithmi¬ 
cally equal interval for r, a linearly equal interval for v. The 
/c^^-bin of the angular momentum slice corresponds to the 
interval [Jma^{k — 1)‘^/Nj , /Nj] and is represented by 

jk = Jmax(fc - 112V INf. 

We modify the splitting algorithm using the fact that 
the angular momentum is an invariant of the Hamiltonian 
system. Hence, one may treat each slice with a different value 
of j in phase-space independently, except for gravitational 
coupling via the Poisson equation. We include the inertial 
component of the force, in the “drift” step (equations]^ 

and|^, while the “kick” step (equation]^ corresponds solely 
to gravitational force: 

r(uuj) = /[r*(-At/2),u*(-At/2),j,t], (9) 

= f{r,v + GMr/r‘^/\t,j), (10) 


/(r,u,j,t + At) = /**[r*(-At/2),u*(-At/2), j], (11) 

where r* and u* solve analytically the motion in absence 
of gravity starting from coordinates (r, u,j) in phase-space 
(see, e.g., Colombi V Toum^|2008 ): 


r*(h) = 


y^2r‘^HK - p + 2 sgn(u)iVK/i 




2Hk 


-u* (A) = sgn(u). / 2Hk - 




(12) 

(13) 


with Hk = v‘^/2 + j^/(2r^) (when u* < 0, these equations 
are valid until v* = 0). 

Because a non-zero angular momentum bends the tra¬ 
jectories in (r, u) space, the drift step requires a two- 
dimensional interpolation of the phase-space distribution 
function in (r, v) space, while the kick step, which only mod¬ 
ifies the velocities, can be completed with a one-dimensional 


interpolation. We follow Fujiwara (1983), and carry out the 


interpolations using third-order splines. In this interpolation 
scheme, however, the positivity of the phase-space distribu¬ 
tion function is not warranted, and numerical aliasing and 
diffusion effects are expected when the phase-space distri¬ 
bution function varies over scales of the order of, or smaller 
than, the mesh element size. 

In order to reduce such numerical artefacts, we modify 
equation as follows: 




Po 


( 27 r cr; 
1 


5 ) 3/2 


exp 


1 u" + fy 




r ^ Rn, 


(14) 


with A = 1/2. Then we recompute po in equation ^ so 
that the total mass remains unity. This apodization slightly 
changes the actual values of the virial ratio to ~ 0.55 and 
0.11, although we shall still denote them by 0.5 and 0.1 just 
for simplicity. It may also modify the long-term dynamical 
properties of the original Henon sphere relative to what is 
expected. This is why we check again the extent to which 
the spherical nature of the system is retained in the A/’-body 
simulations (§[^. 

Adopting a logarithmic binning for r is well suited for 
tracing small-scale features around the center of the system. 
This implies, however, that radii smaller than a finite mini¬ 
mum value Rmin are missing from the computing domain. A 
conventional trick to overcome the problem is to assume a re¬ 
flecting boundary at r = i^min (see, e.g.. Gott 1973 Fujiwara 


1983). Usually, a systematic time-lag between orbits in this 


method is neglected: particles reaching the reflective kernel 
boundary instantly travel the 2Rmin distance through the 
central region, while they should actually take a finite time 
depending on their radial velocity and angular momentum. 
In VlaSolve, we improve the reflecting sphere method by 
taking into account the actual time spent by particles trav¬ 
elling inside the region r ^ i^min, which is made easily pos¬ 
sible by neglecting the gravitational force. Technical details 
about the implementation are provided in Appendix |A1| 

To complete algorithmic details. Appendix [A^ discusses 
the hybrid parallelization of VlaSolve with OpenMP and 
MPI libraries. 
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Nr 

Ny 

Nj 

At 

1024 

1024 

512 

5 X 10-"^ 

512 

512 

512 

10-3 

2048 

2048 

32 

2.5 X 10-"^ 

1024 

1024 

32 

5 X 10-4 


Table 1. The parameters used for the VlaSolve simulations. 


In this paper, we perform 4 simulation runs with differ¬ 
ent resolutions, each for R — 0.1 and 0.5 (Table[^. To cover 
the dynamical range of interest, the computing mesh uses 
Rtciiw — 0.01, -Rmax — 25 aud e-^max — 1.6. T'hc maximum 
amplitude of the velocity is I’max = 2 and 4 for = 0.5 
and 0.1, respectively. With this choice of the parameters, 
the computational domains are sufficiently large to contain 
all the system up to the end of the simulations, which corre¬ 
sponds to t = 100 for R — 0.5 and t = 35 for i? = 0.1. As will 
be illustrated later in phase-space density plots, these final 
epochs are sufficient for the system to have relaxed at the 
coarse level to a meta-stable state through mixing. Strictly 
speaking, this is not the case in the i? = 0.1 case because a 
fraction of the mass escapes from the system (see, e.g.. 


Albada 1982 Sylos Labini 2012 ), as already mentioned in 


the Introduction. 

We adopt a constant time step At throughout each sim¬ 
ulation. Just to stay on the conservative side, we choose a 
resolutely small value of At, despite the increased compu¬ 
tational cost. Note however that excessively small time step 
might artificially increase diffusion effects related to succes¬ 
sive interpolations of the phase-space distribution function 
( Halle|2015 |. 


In Appendix |A3[ a comparison among all the simula¬ 
tions is performed for R — 0.1. It indicates that diffusion 
and aliasing effects discussed earlier are indeed significant, 
despite the apodization of initial conditions, but do not seem 
to affect the dynamical properties of the system. Note that 
is tempting to undersample angular momentum space since 
j is an invariant of the dynamics. However, we show in this 
appendix that it is not wise to do so, because it can provoke 
nonlinear instabilities after a few dynamical times. 


3 iV-BODY SIMULATION WITH GADGET 


We perform the A-body simulations using the latest version 
of the Gadget-2 code ( Springel|2QQ5 ). Only the treecode part 
of this “treePM” algorithm is employed. The particle number 
is varied from N — 10^ to 10^ for = 0.1 and 0.5. We also 
run an additional simulation with A = 10 ® for i? = 0 . 1 . 

We choose the parameters for Gadget runs as follows: 


• The softening length of the gravitational force is set 
as e = 0 . 2 A“^/®, that is about 1/16 of the initial mean 
interparticle distance (47r/3A)^/®i?H (this estimate neglects 
the effects of the apodization . 

• In Gadget, each particle has its individual time step 
bounded by dt — min[dtmax, ( 2 ? 7 e/|a|)^/^], where a is the 
acceleration of the particle and 77 is a control parameter. We 
choose ?7 = 0.025 and Atmax = 0.01. 

• The tolerance parameter controlling the accuracy of 
the relative cell-opening criterion (parameter designed by 


ErrTolForceAcc in the documentation of Gadget, see equa¬ 
tion 18 of Springel||2005 ) is set as ap = 0.005. 

Appendix presents the effects of changing these pa¬ 
rameters on the phase-space distribution function for simu¬ 
lations with A = 10® particles and a virial ratio of = 0.1. 
These analyses, performed at t = 15, confirm that the pa¬ 
rameters used for the simulations of this paper are reason¬ 
able. Interestingly, changing the softening length by large 
factors does not influence much the results, as already no¬ 


ticed previously in the literature (see, e.g. Barnes, Lanzel, 
fc Wiffiam^|2009 ), as long as it is kept small enough. 


4 CONSISTENCY CHECK: SPHERICITY OF 
THE A-BODY RESULTS 


Before presenting comparisons between Gadget and 
VlaSolve, it is necessary to make sure that the sphericity 
of the system is preserved in the Gadget simulations be¬ 
cause our Vlasov runs are performed assuming exact spher¬ 
ical symmetry. Figure shows, for different values of the 
number of particles A, the evolution with time of the ratios 
b/a and b/c, where a ^ b ^ c are the eigenvalues of the 
inertia tensor of the particle distribution. 

The dashed regions correspond to the one sigma zone 
obtained from an ensemble of 100 local Poisson realizations 
of the spherical density p(r), which is estimated from in¬ 
terpolation over spherical shells from the Gadget particles. 
From the measurements in Fig. deviations from spheri¬ 
cal symmetry due to the particle shot noise can be roughly 
scaled to 


where and are the variances of 6 /a and b/c ob¬ 
tained from the dispersion over the 100 realizations. Note 
that equation (15) is not intended to be accurate. The as- 
phericity due to discreteness should depend on details of 
the density profile, as shown in Fig. While it would be 
possible to compute in a perturbative way the quantities in 
equation (15) from statistical analysis of the inertia tensor 
assuming N ^ 1 and using error propagation formulae, this 
is a cumbersome exercise far beyond the scope of this paper. 

We also note that another possible source of errors 
comes from the position of the center of the system. In¬ 
deed, an inaccurate determination of the center obviously 
worsens the apparent agreement with spherical symmetry. 
In the measurements presented in Fig.[^ the inertia matrix 
is not computed with respect to the center of gravity of the 
particle distribution, which can be affected by the fact that 
some particles can get far away from the system through 
A-body relaxation. Instead, we determine the center of the 
system using an iterative procedure trying to optimize the 
match of the phase-space distribution function with that of 
the Vlasov code, as detailed in § |6.1| This procedure is not 
free from errors either, and may contribute to the fluctua¬ 
tions observed in the curves of Fig. 

Inspection of Fig. shows that the measured ratios b/a 
and b/c behave differently in the R = 0.5 and R = 0.1 simu¬ 
lations. In the R = 0.5 case, the agreement of the measure¬ 
ments with the Poisson prediction is in general good, with 
a slight trend to ellipticity, except for the top red curve and 
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Figure 1. Evolution of the departure from spherical symmetry: ratios of the eigenvalues of the inertia tensor of the system as functions 
of time in the Gadget simulations. To emphasize the small differences from unity, the quantity sgn(l — Taxis) logio 1^ “ ^axisl is plotted 
as a function of time, where Taxis = bja (upper curves on each panel) or b/c (lower curves) and a ^ b ^ c are the eigenvalues of the 
inertia tensor of the Gadget particle distribution. Each color corresponds to a given value of the number N of particles as indicated in the 
panels. Dashed regions correspond to the one sigma confidence level zone expected for a particle distribution locally Poisson sampling the 
spherically symmetrical projected density profile p(t, t), where p{r,t) is estimated from interpolation of the Gadget particle distribution 
in spherical shells. To calculate the average of Taxis cind the associated one sigma error contours, 100 local Poisson realizations have been 
performed for each snapshot and value of N considered, except for N = 10^ and N = 10® (on right panel only for the latter). In the last 
cases, the dashed regions correspond to an extrapolation of the results obtained from N = 10®. 


the bottom green curve where the deviation from spherical 
symmetry is larger than the Poisson expectation. Still, in 
the case of R = 0.5, the system remains to a very good ap¬ 
proximation spherical for all values of N, given the expected 
deviations due to pure statistical noise. 

The curves representing the eigenvalue ratios are more 
steady for = 0.1 than for R = 0.5, which might be slightly 
puzzling at first sight. However, a very plausible explana¬ 
tion of this difference is that the initial velocity dispersion 
is larger for R = 0.5 than for = 0.1, hence adding a more 
prominent random component to the time behavior of the 
deviation from sphericity. 

Regarding R = 0.1, deviations from spherical symme¬ 
try are clearly more significant compared to local Poisson 
expectations after t ^ 3, roughly the collapse time of the 
sphere. While the iV = 10^ run exhibits a deviation larger 
than 10 percent, spherical symmetry is confirmed to be a 
good approximation for N ^ 10^. 

Finally, we also check deviations from spherical sym¬ 
metry for subsets of particles in excursions corresponding 
fo / ^ /th, where / is the phase-space distribution func¬ 
tion measured in the 1024 x 1024 x 512 VlaSolve runs. For 
each value of the virial ratio, two thresholds /th are chosen 
such that the excursions contained initially about 90 and 
60 percent of the total mass (see bottom panels of Fig. 
below). Given the uncertainties in the measurements, the 
above conclusions still hold: the properties of the deviations 
from spherical symmetry, that we do not show here for sim¬ 
plicity, do not indeed depend significantly on radius. We 
only notice a slight improvement in the R = 0.5 case when 
considering particles in the excursions. 


5 PHASE-SPACE DENSITY: VISUAL 

INSPECTION 

Now we are ready to perform direct comparisons between the 
Vlasov and iV-body simulation results. For this purpose, we 
consider the phase-space density at different epochs (Figs. 
tobelow). To be more specific, we plot the constant angu¬ 
lar momentum slice of f{r,v,j) at j = 0.244, and its integral 
over the angular momentum: 

fsummedVv) = J f {r, V, j) 2njdj. (16) 

Figures[2]and|^plot f{r,v,j ~ 0.244) and /summed(r, t>), 
respectively, for the VlaSolve and Gadget simulations of the 
warm case, R = 0.5. In both figures, snapshots at t = 10, 
50, 80 and 100 are plotted from left to right. The pan¬ 
els correspond to the VlaSolve runs with {Nr, Ny, Nj) = 
(2048, 2048, 32) and (1024,1024, 512), the Gadget runs with 
N = 10^, 10® and 10®, from top to bottom. 

The overall conclusion of the visual inspection of Figs.|^ 
and is that the Vlasov solver and the A-body code exhibit 
very good agreement with each other, probably even much 
more than expected. In particular, both results present a re¬ 
markably similar instability in the region 1 t 10®‘®, even 
in details, showing a surprising reliability of the conventional 
A-body approach for these particular initial conditions. 

However, before reaching this conclusion, one has to 
take into account several limiting factors. In particular, we 
should bear in mind the fact that the VlaSolve simulations 
are subject to significant diffusion, which smears out fine de¬ 
tails of the phase-space distribution function. This diffusion 
effect is clearly visible at t = 50, when comparing the outer 
filamentary structures observed in the Vlasov simulations to 
the A-body result. Putting aside this coarse-graining effect, 
the structures are exactly similar in both the A-body and 
Vlasov simulations at t ^ 50, even including small gaps in 
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Figure 2. VlaSolve versus Gadget in phase-space: phase-space density for R = 0.5 and averaged over j £ Ij = [0.225,0.264]. Each 
column of panels corresponds to a given value of time t, increasing from left to right. The first two lines of panels display f{r,v,j) for 
VlaSolve simulations with (Nr, Ny, Nj) = (2048, 2048, 32) and (1024,1024, 512) respectively, while the three bottom lines correspond to 
the A/’-body simulations, with various values of the number of particles N as indicated on each panel. Note that the VlaSolve simulation 
with (Nr,Nv,Nj) = (2048,2048,32) has only one angular momentum slice, J = 0.244, in the interval Ij, so there is no blurring of 
the filamentary details of f(r,v,j) on the left side of the peak of the distribution function contrarily to the other cases. In the A^-body 
case, /(r, v,j) was computed on the same mesh as the (1024,1024, 32) VlaSolve simulation using nearest grid point interpolation, which 
explains the artefacts on the color pattern in the last two lines of panels. ^ 2014 RAS MNRAS 000 p~[fT^ 
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Figure 3. Same as in Fig.|^ but the phase-space distribution function has now been summed up over the whole available range of values 
of j G [0, Jmax = 1-6], where Jmax is the maximum sampled value of j for the VlaSolve simulations. 
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Figure 4. Same as in Fig. but for a colder initial configuration with virial ratio R = 0.1. There is also an additional line of panels 
corresponding to the Gadget simulation with N = 10® particles. Note the large R tail escaping from the system, corresponding to a 
fraction of the mass with positive energy (see, e.g., |van Albada| 1982 1 [Joyce, Marcos, fc Sylos Labini|200^ |Sylos Labini|2012| >. 


the phase-space distribution function related to nonlinear 
instabilities that start building up. These instabilities grow 
further at later epochs. They are considerably smeared out 
in the (1024,1024, 512) VlaSolve simulation but unquestion¬ 
ably present. Adding resolution in (r, v) space (at the cost of 
resolution in j) improves the agreement with Gadget, which 
confirms that the instabilities observed in the Gadget simu¬ 
lations are physical and not of numerical nature. 

Figure [^indicates that lowering the number of particles 
in the A^-body simulations may be interpreted as a coarse- 


graining: it makes finer details more fuzzy but still keeps 
global features of the phase-space density correctly. We also 
note that using a small number of slices in j in the Vlasov 
solver does not seem to alter the dynamical properties of 
the system despite the considerable level of aliasing it intro¬ 
duces. 

The situation is more complicated for the cold case, 
R = 0.1 (Figs. and 1^. Up to t gy 10, the above conclu¬ 
sions for R = 0.5 are still valid. However, some instabilities 
emerge at t ~ 10 in the Gadget simulations with V ^ 10® 
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Figure 5. Same as in Fig. but for R = 0.1 and with the additional A/’-body simulation involving N = 10® particles. 


particles as well as the (2048,2048,32) Vlasov run. Until 
this epoch, the N ^ 10^ and the (1024,1024,512) simula¬ 
tions agree perfectly with each other (modulo the smearing 
effects already discussed above) and present a smooth phase- 
space density without any sign of instability. On the other 
hand, the other simulations exhibit slightly irregular phase- 
space density. Such a trend is easily seen in Fig. even 
though not so obvious in the (2048,2048,32) Vlasov simu¬ 
lation. These irregularities appear as well in the N ^ 10^ 
simulations but at later epochs, and then develop in a dra¬ 
matic way. A careful inspection of successive snapshots of 


the simulations indeed suggests that the onset of these ir¬ 
regular patterns comes later with increasing N. 

As discussed in Appendices [B]and |A3| these instabilities 
result from the discrete nature of the system in the V-body 
case, and from the aliasing effect due to sparse-sampling of 
the angular momentum space in the Vlasov code. Since the 
pattern of the instabilities changes significantly from one 
simulation to another unlike the i? = 0.5 case, they should 
be due to numerical, not physical, origin. 

As shown in Appendix their presence is very insen¬ 
sitive to the choice of softening, time step or parameters 
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controlling force accuracy in Gadget. They can therefore be 
reduced only by increasing the number of particles and the 
resolution in the Gadget and VlaSolve simulations, respec¬ 
tively. 

It is important to notice that even the N = 10^ re¬ 
sult might be insufficient to describe properly the system at 
late epochs. In the (1024,1024, 512) Vlasov simulation, the 
phase-space distribution function seems to be rather smooth 
at all times and the system is free of instability, contrarily 
to the = 0.5 case. However, it is difficult at this point to 
know if actual physical instabilities build up at late times in 
the = 0.1 case, because diffusion in the Vlasov simulation 
might prevent the appearance of some unstable modes. 

While the irregular patterns observed in Figs. and 
are definitely of numerical nature, the fact that they develop 
so easily may indicate that the system is prone to react non¬ 
linear ly to small perturbations. Uneven gaps between the fil¬ 
aments of the phase-space density can be observed at t = 15 
(third column of Fig.[^, even in the (1024,1024, 512) Vlasov 
simulation, and one might expect that they correspond to 
seeds of actual physical instabilities. In this respect, the sys¬ 
tem might actually develop, at some point, physical unstable 
modes. These results are quite suggestive of what was ob¬ 
tained previously with a spherical shell code for cold and 
self-similar systems ( Henriksen V Widrow||1997 ). 

Even with our V = 10^ particle simulation, it is not 
clear whether these unstable modes dominate over collective 
effects due to discreteness. A better understanding of the 
phenomenon would require a convergence study using even 
higher-resolution simulations. 


6 STATISTICAL ANALYSIS 

6.1 Correlators and entropic estimators: 
definitions and concepts 

To perform a more accurate analysis, one can try to quantify 
to which extent the particle distribution in the V-body sim¬ 
ulations can be considered as a local Poisson process of the 
phase-space density calculated in the semi-Lagrangian code. 
To do so, we use, in addition to entropic measurements de¬ 
scribed further, the following correlators, 

Ck = —, (17) 

with 

N 

i=l 

Kk = (19) 

In these equations, /c is a positive integer, / the VlaSolve 
phase-space density, M the total mass, dQ = 27rdr x dvxjdj 
and = (n,'^’i,A), where Pi^ Vi and ji are respectively the 
radial position, radial velocity and angular momentum of 
each particle of the Gadget simulation. 

For a point set randomly sampling a smooth density 
distribution g, the probability density p(f^) of having a given 
particle at phase-space position Q is independent from the 
rest of the particle distribution and is simply proportional 
to s(0): 


p{n)dn = (20) 

The density probability of having N particles at respective 
positions Ui, U 2 , ..., 0.n is given by 

N 

p(ni,...,njv) = Hp(^»)- (21) 

i=l 

Ensemble averaging of gk under the law g then reads 


{f^k)g — 


n N 


7^(Qi,...,Qiv) dQi---dUiv, 


M 

TV 


r 1 


—[/(n,)]''s(n,)dn,, 


n/ 


M 


-dQn 


-I 


[/(n)]'=5(n)dn, 


and 


{Ck)g — 


f[fiQ)fg{Q)dQ 

/[/(n)]'=+idn ■ 


( 22 ) 

(23) 

(24) 

(25) 


Hence, if the distributions g and / coincide, i.e., in our case, 
if Gadget actually Poisson samples the VlaSolve phase-space 
density, one obtains {Ck)g=f = 1 after ensemble averaging. 

When increasing /c, more weight is given to regions in 
phase-space corresponding to larger values of /. Eor a point 
process totally anticorrelated with /, Ck cancels, while its 
largest possible value is given by Ck = (Mmax/^)//^/c > 1, 
when all the particles stay in the region where / is maximal. 

An important issue is to compute properly the center of 
the system position in the Gadget simulations. In order to 
do this, we find the coordinate origin maximizing Ci, even 
though the result of such a procedure can potentially lead 
to Cl > 1, to optimize the match between concentrations of 
particles and local extrema of /. 

The variance of Ck can also be calculated in an analo¬ 
gous way to {iik)g- 


^Cl = {Cl),-{Ck)l 


Af / V 1 , . 2 

-J^{g2k)g - -^{Wg 


N 


(26) 

(27) 


which reduces to AC^ = (M/A)(/^ 2 fc//^fc) — 1/A when / and 
g coincide. In practice, we shall use the following estimator 
for this statistical error: 


AC," 


'M 1 2 ' 

iV ? 


(28) 


where 112 k and /x, are directly estimated from the A-body 
simulation. 

The log-likelihood that the Gadget particle distribution 
locally Poisson samples the VlaSolve phase-space density / 
can be written, following the reasoning that leads to equa¬ 
tion (21), 


In C = In 


fm 


M 


(29) 
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However, the region D where / > 0 being of finite extent, one 
expects ln£ = — oo as soon as a particle escapes P, which is 
very likely, due for instance to A/’-body relaxation. Further¬ 
more, the Vlasov solver does not guaranty the positivity of 
/. To take into account in a fair way both the defects of the 
V-body and the Vlasov simulations, it is better to restrict 
to a region Pth where / is strictly positive: 


Dth = such that f{fl) > fth, fth > 0}. (30) 


The log-likelihood of having Q ^ N particles in the region 
Pth and the rest outside it (leaving the freedom of the re¬ 
maining particles to span all the space outside Pth) is given 
by a binomial law: 


ln£b(Q, = In 


V! 




(1-^) 


N-Q 


(31) 


_{N-Qy.Q\ 

where ly is the fractional mass inside Pth in the VlaSolve 
simulation. Hence, equation (29) simply becomes 

'fm' 


In £ = In 


Mth 


+ In ChiQth, v), 


(32) 


between the measured value of the log-likelihood and its ex¬ 
pectation under the law /: 


5S= — [{ln£}f-\nC], 


(38) 


where C is given by expression (32) calculated for Qi ex¬ 


tracted from a Gadget simulation. The quantity 6S esti¬ 
mates the magnitude of the difference between the under¬ 
lying smooth phase-space density g sampled by Gadget and 
the VlaSolve phase space density, /: its ensemble average 
other many Gadget realizations indeed reads, when neglect¬ 
ing the binomial term in equation (|32|), 


(SS), ^ 


/ 


^[5(n)-/(n)]in 


m 

Mth 


dft. 


(39) 


Under the assumption that the V-body simulation Poisson 
samples the distribution /, the magnitude of (5S' should be 
of the same order of cr^. 


6.2 Correlators and entropic estimators: 
measurements 


where Qth is the number of particles of the Gadget simula¬ 
tion inside Pth and Mth = dP f(P). 

Note that the distribution of particles which maximizes 
the first term in equation ( [^ corresponds again to the case 
where all the particles of Pth stay in the region where / 
is maximal, similarly to the case when the correlator Ck is 
equal to its maximum possible value. Clearly, this situation 
is not typical, but it is in fact the most likely to consider 
when it can take place: this is why we maximize Ci to esti¬ 
mate the center of the V-body system, even though it might 
turn to be larger than unity. 

The expectation value of ln£ under the law / can be 
obtained by ensemble averaging: 

5(/th) = -^(ln£>/ = 5/(/th) + 5b(/th), (33) 

N 

Shifth) = (35) 

Q=0 


m) 


Mth 


dP, 


(34) 


In the limit /th ^ 0, the quantity Sf{fth) reduces to the 
Gibbs entropy of the system, which explains the choice of 
notations. Moreover, if V ^ 1 and if the fractional mass u 
inside the domain of interest Pth is of order of unity, which 
is the case for our analyses, the term *S'b(/th) is in prac¬ 
tice negligible compared to 5'/(/th), so S{fth) depends only 
weakly on the total number of particles, as expected. 

The variance of In C can be calculated likewise 


2 


ctl 


yVb [(ln£^)/-(ln£)^] 


1 (/ /(«),„, 

\fm] 

IP, Mth 

Mth 


(36) 

dn-i.[5(/th)]'|, 

(37) 


where we have neglected, following the arguments developed 
earlier, the contributions of Sh to the error. 

To understand better the interest of using the statis¬ 
tics given by equation (32), one can introduce the difference 


Top panels of Fig.|^show the quantity Sf{fth) as a function 
of time for the various VlaSolve simulations we performed 
and two values of fth chosen such that approximately 90 
percent and 60 percent of the total mass is initially inside 
the excursion Pth, respectively. The quantity 5'/(/th) is a 
Casimir invariant -that is an integral over a function of /- 
and should thus be conserved during runtime if the code was 
perfect. This not the case because of diffusion and aliasing 
effects in (r, v) space: deviation from conservation of Sf hap¬ 
pens shortly after collapse time. Then there is a strong mix¬ 
ing phase during which Sf increases, then possibly decreases, 
according to the value of /th, and finally reaches an approx¬ 
imate plateau. Deviation from conservation of Sf naturally 
happens sooner when resolution in {r,v) space is smaller. 
Resolution in j space does not have much influence on Sf 
because angular momentum is an invariant of the dynamics. 
However, as clearly shown in § [^ and in Appendix |A3| for 
P = 0.1, we already know that sparse sampling in j space is 
not recommended since it can introduce some instabilities in 
the dynamics, even though this effect does not affect much 
our likelihood measurements. 

Middle panels of Fig. show the quantity — ln£/ (Nu) 
measured in Gadget from the particles belonging to the ex¬ 
cursion Pth as a function of time, where ln£ is given by 
equation (32). For a given value of the threshold /th, if the 
Gadget simulations would actually behave like Poisson real¬ 
izations of the VlaSolve ones, all the colored curves should 
be close to the solid line, which corresponds to Sf. This is 
clearly not the case for small fth (upper group of curves), ex¬ 
cept a early times. Increasing the number of particles in the 
iV-body simulation improves the agreement with the Vlasov 
code for R = 0.5 but does not seem to have a convincing 
impact in the P = 0.1 case: for fth = 0.02, all the N- 
body simulations converge to the same plateau somewhat 
below the Vlasov code result. On the contrary, for fth = 0.2 
and P = 0.1, the agreement between Gadget and VlaSolve 
is striking at all times, except may be for the V = 10^ 
simulation during the strong mixing phase. Note also, that 
at late times, all the V-body simulations converge which 
each other, independently of fth and P, except again for the 
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Vlasov entropy, / ^0.0025 k. 0.02, R=0,5 


Vlasov entropy, / ^0.02 Sc 0.2, /? = 0.1 






Vlasov vs 7V-body, / ^0.0025 Sc 0.02, P=0.5 Vlasov vs TV-body, / ^0.02 &: 0.2, /?=0.1 




Figure 6. Entropic measurements: effects of VlaSolve resolution (top two panels) and Gadget number of particles (four bottom panels). 
The left and right panels correspond respectively to R = 0.5 and 0.1. On the top panels, the quantity S'/(/th) given by equation ( |34| ) is 
plotted as a function of time for the Vlasov simulations and for two values of /th indicated on each panel corresponding to approximately 
initially keeping 90 and 60 percent of the mass inside the excursion. Each curve corresponds to a given resolution as indicated on each 
panel (the dashes are nearly superposed to the solid line). The top/bottom group of four curves correspond to a smaller/larger value 
of /th- On the middle panels, the solid line is the same as on the top panels, while the colored curves display, for each value of the 
particle number N in the Gadget simulations, the quantity —\n/l/{Nh') as a function of time, where \nC is given by equation ( |32| ). 
If the V-body simulations would Poisson sample the VlaSolve phase-space density, the ensemble average of this quantity over many 
Gadget realizations should match the solid line (except for a negligible correction due to the S\^ term in equation |34| ). Finally, the bottom 
panels show the fractional mass as a function of time for the two values of /th considered. On the two bottom right panels, there is an 
additional purple curve nearly indistinguishable from the red one, corresponding to the additional simulation with 100 millions particles 
we performed for = 0.1. In the four bottom panels, the thickness of each colored curve takes into account statistical errors (equation 
|37| for ln£). In addition, for the middle panels, systematic errors due to the interpolation of the phase-space distribution function in the 
VlaSolve simulations also contribute to the estimated errors. In the latter case, we compute f(^i) both using nearest grid point and 
linear interpolation from the values of / on the computational mesh. The difference between the two interpolating methods adds to the 
thickness of the curves. Note that we use the (1024,1024, 512) VlaSolve simulation to perform the comparison to V-body results, to 
minimize the effects of interpolation. 
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N = 10^ simulation with R = 0.1, but we know that this lat- 
ter presents significant deviations from spherical symmetry 
and should be probably discarded for the analyses performed 
here. 

To complete the analyses and understand better the 
results obtained for the log-likelihood, the fractional mass 
inside the excursions / ^ /th is shown in bottom panels of 
Fig. Again, this quantity is a Casimir, so it should not 
change with time in the idealistic case. In practice, while it 
is difficult to predict the effects of aliasing on the VlaSolve 
mass inside Dth, diffusion effects are more likely to decrease 
it, especially by dilution of filamentary structures that build 
up during the course of dynamics. In the R = 0.5 case, 
most of the disagreement between Gadget likelihood and its 
expectation given by VlaSolve can be understood in terms 
of fractional mass: effects related to the discrete nature of 
the A’-body simulations seem to spread particles away from 
Hth- However this process is subtle and seems to remain local 
as suggested by visual inspection of Figs. and We also 
checked that it does not affect dramatically the projected 
density, p{r). 

In the = 0.1 case, the interpretation of the results is 
slightly more complicated. For /th = 0.2, the Gadget frac¬ 
tional mass inside the excursion Dth behaves similarly as 
in the i? = 0.5 case as a function of particle number. On 
the other hand, when examining the quantity — ln£/(Az/), 
the A-body measurements converge with each other and 
with VlaSolve much better, especially after relaxation. This 
means that particles left in Dth are redistributed in a non 
trivial way, such that the effects of the excursion mass loss 
are compensated. For /th = 0.02, even the A = 10® Gadget 
sample disagrees with the VlaSolve simulation. Clearly, 
the Vlasov simulation becomes quickly defective in regions 
where / is small. On the other hand, convergence of the 
Gadget simulations at late times might be misleading. In¬ 
deed, we noticed from visual inspection of Figs. and 
that some instabilities appeared in all of them as soon as 
t ^ 15, although later when A is larger. Interestingly, the 
measurements in the A = 10^ and A = 10® simulations 
are nearly indistinguishable from each other, which is a sign 
that we are nevertheless close to numerical convergence. 

Entropic measurements of Fig.j^are confirmed, at least 
partly, by Fig. HI In particular, a depression of which the 
depth depends on the number of particles in the A-body 
simulation appears on all the curves. When increasing A, 
the amplitude of the depression decreases and the occur¬ 
rence of its maximum amplitude is delayed, independently 
of the actual dynamical state of the system. Again, it can 
certainly be attributed to collective effects due to Poisson 
noise. Overall agreement between A-body and Vlasov codes 
improves when increasing the number of particles in the A- 
body simulation. For = 0.5, this is rather independent of 
k in equation (17), i.e. of the fact of putting more or less 
weight to overdense regions in phase-space. In the i? = 0.1 
case, putting aside the depression of which the depth de¬ 
pends on the number of particles, the correlator Ci starts to 
decrease with time at t ~ 10. This can be mainly attributed 
to defects in the Vlasov simulation in underdense regions 
as discussed earlier. For k ^ 2, which gives more weight to 
higher values of the phase-space density, the correlator in¬ 
deed stays steady as a function of time (again putting aside 
the A-dependent depression). However, one notices for /c = 3 


a net increase with time of the correlator for the simulation 
with A = 10^ particles, but let us remind that this simula¬ 
tion presents significant deviations from spherical symmetry. 


7 CONCLUSION 


In this paper we have compared the phase-space distribution 
function traced by the particle distribution in Gadget sim¬ 
ulations to the results obtained with our new Vlasov code 
VlaSolve for spherical systems, an improved version of the 
splitting algorithm of Fujiwara] ( |1983| . For the specific com¬ 
parison, we have chosen (apodized) Henon spheres, which 
are known to be insensitive to radial orbit instability and 
in particular to preserve the spherical nature of the system. 
The latter property is confirmed from simulations run with 
three-dimensional A-body codes. We considered two values 
of the initial virial ratio of the spheres, = 0.5 and = 0.1, 
corresponding to “warm” and “cold” configurations, respec¬ 
tively. 

We have plotted detailed structures of the phase-space 
distribution functions varying the spatial/mass resolution 
of the numerical code in a systematic fashion, we have con¬ 
ducted further a quantitative analysis by introducing two 
new statistical tools. The first one is of entropic nature and 
corresponds to the log-likelihood quantifying to which extent 
the A-body results represent a local Poisson sampling of the 
Vlasov phase-space density. The second tool is a correlator 
of order k, proportional to the integral over phase-space of 
the product between the Vlasov phase-space density raised 
to the power k and the particle distribution function. 

The overall conclusion is that both the Vlasov and A- 
body methods agree remarkably well with each other, both 
from the visual and statistical points of view, if sufficient res¬ 
olution is employed. Given the completely different numer¬ 
ical approaches to collisionless dynamics, this is not trivial 
at all, and the degree of agreement that we have shown for 
the first time is perhaps even better than what had been 
expected before. This is reassuring for numerous previous 
results that have been almost exclusively obtained from the 
A-body method. 

Nevertheless there are still unsolved subtle issues in de¬ 
tails: 


• When performing a visual inspection of the phase-space 
distribution function in the cold case, = 0.1, although still 
good at the coarse level, we find that the level of agreement 
between the A-body and the Vlasov codes worsens at small 
scales after a few dynamical times. This is mainly due to col¬ 
lective effects induced by the shot noise of the particles in the 
A-body simulations (and not to close particle encounters). 
Even with A = 10® particles, we are not able to prove nu¬ 
merical convergence of the A-body results. The comparison 
at this level, however, is made difficult by the fact that the 
Vlasov code is significantly diffusive, which might prevent 
the development of a variety of physical unstable modes. 

• While the statistical tools do not provide as rich and in¬ 
tuitive information as visual inspection, they identify some 
subtle effects. In particular, when taking into account gen¬ 
eral trends due to diffusion in the Vlasov code, significant 
for R = 0.1, we notice that the match between Gadget and 
VlaSolve worsens with time, then improves. The degree of 
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Figure 7. Correlators between VlaSolve and Gadget as functions of time. These quantities, defined in equations ( |17| ), ( |18| ) and ( |19| ), 
are plotted for = 1, 2, 3 increasing from top to bottom, while left and right panels correspond to = 0.5 and 0.1, respectively. The 
thickness of the curves, analogously to Fig. takes into account statistical errors according to equation ( |28| ) using the measured value 
of v'2k cind Uk cind systematic errors due to the interpolation of the phase-space density in the VlaSolve samples. Note that there is an 
additional purple curve on each panel of the right column corresponding to the 100 millions particles simulation. 


the mismatch increases, and it shows up earlier, when re¬ 
ducing the number of particles in the A/’-body simulation. 
Again, this may be ascribed to collective effects due to the 
shot noise of the particles. Nevertheless, the very good match 
between the Gadget simulations with N — lY and N — lY" 
particles may suggest that convergence is nearly reached in 
terms of number of particles and information theory, even if 
it is not fully proved. 


It is worth mentioning again that the collective effect 
mentioned above is not related to A/’-body relaxation, but 
rather results from random Poisson fluctuations. This can 
be formulated as follows (se e [Aarseth, Lin, fc Papaloizou 


1988| [Henriksen Widrow| 1997[ |Boily, Athanassoula, fc 
Kroupa|2002[ [Joyce, Marcos, fc Sylos Labini|200^ for simi¬ 

lar arguments): a given particle at some distance r from the 
center of the system feels a force proportional to the num- 
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ber A^in of particles inside the sphere of radius r. Poisson 
fluctuations imply thus that there is a relative error of order 
of 1/y/Nin on this force. Importantly, the inner number of 
particles Nin changes with time with random fluctuations 
around the mean behavior: these fluctuations can be con¬ 
sidered as a correlated random walk. Indeed, because of the 
finite velocity dispersion, particles cross both inwards and 
outwards the frontier of the sphere of radius r. A larger ve¬ 
locity dispersion weakens the amount of correlation, thus 
makes the errors on the force more random, which should 
have a fuzzy effect on the phase-space density, similarly as 
collisional relaxation: this is what we can expect for R = 0.5 
and as observed on Fig.[^ On the contrary, a smaller veloc¬ 
ity dispersion makes the error on the force more systematic 
which should induce coherent distortions of the phase-space 
density: this is what we can expect for = 0.1 and con¬ 
firmed by visual inspection of Fig. This effect has non¬ 
trivial consequences on the energy spectrum of the particles, 
particularly in cold configurations ( [Joyce, Marcos, Sylos 


Labini|20Q9 ). It certainly explains as well the deviations be¬ 


tween VlaSolve and Gadget observed when measuring the 
statistical estimators defined in this paper. According to 


Aarseth, Lin, & Papaloizou (1988), this collective effect is 


dominant over A’-body relaxation, and, as confirmed by our 
detailed numerical tests in Appendix is not significantly 
influenced by softening. 

Note as well that shot noise creates anisotropies in the 
system, i.e. deviations from spherical symmetry that may 
be eventually amplified. Aarseth, Lin, & Papaloizou (1988) 
argue that this effect is subdominant compared to the radial 
component of the noise-induced perturbation when consider¬ 
ing the collapse of an homogeneous sphere. Although their 
calculation is performed only prior to collapse and in the 
cold case, we believe that the conclusion still remains valid 
for the kind of initial conditions studied in this paper, as sug¬ 
gested by our numerical experiments that seem to preserve 
well spherical symmetry. 

Clearly, the collective effect due to particle shot noise 
is a real problem for simulations of close to cold spherical 
systems when it comes to examine fine structures of the 
phase-space density. We were not able to prove convergence 
of the phase-space density in the R = 0.1 case even for 
an A = 10^ particle simulation. Notably, this may have 
non-trivial consequences on the fine structure of simulated 
dark matter halos, where numerical convergence in terms of 
number of particles might not have been reached yet despite 
the numerous intensive studies. Indeed, convergence toward 
the continuous limit might be much slower than expected, 
hence giving the false impression that it is achieved. 
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most expensive part of the code, can thus be computed in¬ 
dependently for each slice of constant j. We therefore easily 
reach an almost perfect parallelization up to a number of 
tasks equal to the grid resolution Nj of angular momentum 
space, which is typically larger than the number of available 
cores on a shared memory system. 

Distributed memory parallelization via MPI is not as 
simple. Indeed, spline interpolations are intrinsically non¬ 
local, which makes the parallelization along dimensions 
other than j non trivial. Sticking with the trivial paralleliza¬ 
tion described above unfortunately limits the maximum to¬ 
tal number of processes running in parallel to Nj, which 
is suboptimal. We overcome this limitation by performing 
MPI domain decomposition in (r, v) space, following the ap¬ 
proach of Grouseilles & al. (2009), who propose to localize 
the cubic spline interpolation to each domain by using Her- 
mite boundary conditions between the domains with an ad 
hoc reconstruction of the derivatives. 


APPENDIX A: VLASOV SOLVER: DETAILS ON 
THE ALGORITHM 

Al Reflecting boundaries with time delay 

In this appendix, we explain how reflecting boundaries con¬ 
ditions with time delay are implemented in VlaSolve. 

If the mass inside the sphere of radius Rmin is neglected, 
the trajectories followed by each test particle associated to a 
grid site that penetrates the sphere are fixed and do not de¬ 
pend on time. This property, combined with the fact that we 
use a constant time step, allows us to pre-compute these tra¬ 
jectories once and for all. The delayed central sphere method 
is then implemented by associating a linked list to each grid 
site whose associate test particle radial position r half a time 
step backward in time is such that r ^ Rmin- Each linked 
list contains as many elements as the number of time steps 
needed for the particle to travel a distance of 2Rmin and 
the element in the list stores the coordinates of the test 
particle n time steps backward in time. Before starting the 
simulation, we initialize each element coordinate and the 
corresponding value of the initial distribution function. For 
each time step, the value of each element is then simply up¬ 
dated by assigning to it the value of its successor while the 
last element value, whose coordinates fall inside the com¬ 
puting domain, r ^ Rmin, is interpolated. A comparison of 
the results obtained with the reflective central sphere to our 
improved delayed central sphere is shown on figure [Al] The 
improvements are unquestionable. 


A2 Parallelization issues 

We implemented a hybrid shared and distributed memory 
version of VlaSolve via the OpenMP and MPI libraries, 
respectively. 

Shared memory parallelism is relatively straightforward 
to achieve in the spherically symmetric case, by taking ad¬ 
vantage of the fact that the angular momentum j is a con¬ 
served quantity. Spline interpolations, which represent the 


A3 Effects of resolution 

Figures |A2| and |A3| show, respectively for j = 0.244 and 
integrated over angular momentum, the phase-space distri¬ 
bution function measured in VlaSolve simulations with dif¬ 
ferent resolutions. These simulations have been performed 
for a Henon sphere with initial virial ratio R = 0.1. Beside 
the very good global agreement between the various runs, 
these figures bring out three effects, which increase when the 
resolution of the phase-space grid is reduced: 


• Diffusion smearing out fine details that build up in 
phase-space during the course of dynamics, for instance 
clearly visible when one compares top to bottom middle pan¬ 
els of Fig. |A2| One concern with diffusion is that it might 
prevent the appearance of unstable modes. However, we did 
not perform any simulation in this work that would prove 
this. 


• Aliasing due to artificial oscillations in the spline in¬ 
terpolation: for the problem studied here, aliasing becomes 
particularly visible after relaxation in the region above the 
large r tail, but this does not have significant impact on the 
dynamics. 

• Aliasing due to undersampling angular momentum 
space: it is visible at all times when one examines the phase- 
space distribution function integrated over angular momen¬ 
tum (top panels of Fig. A3) and can have dramatic con¬ 
sequences on the dynamics. The two top lines of panels of 
Fig. |A2| and |A3| corresponding to a sparse sampling in j 
space with only 32 slices, indeed show the appearance of 
an instability, which presents, on the third column of these 
figures, the same pattern whether (Ar,W) = (2048,2048) 
or (1024,1024). This instability is not present in the sim¬ 
ulations with higher resolution in j, as shown by the two 
bottom lines of panels. Note that the presence of this insta¬ 
bility depends on initial conditions: for R = 0.5, we did not 
notice it for the time coverage considered, t ^ 100 (upper 
line of panels of Figs. and [^. 
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Figure Al. Comparison between the reflecting central sphere method (left panel) and our improved delayed central sphere implemen¬ 
tation (right panel). A simulation of a Henon sphere with (Nr, Ny, Nj) = (200, 200, 200) and a virial ratio R = 0.5 is shown at t = 30 in 
the (r,u = 0,j) plane. The systematic artificial speed increase undergone by orbits that penetrate the central region compared to their 
higher angular momentum counterparts can clearly be observed at low j on the left panel where a reflective sphere is used, while the 
distribution function does not exhibit such spurious features when a delayed kernel is used (right panel). 
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Figure A2. Effect of resolution in the Vlasov code: phase-space density for = 0.1 and j = 0.244. Each column of panels corresponds 
to a given value of time t, increasing from left to right, while each line correspond to a given resolution, (Nr, Ny, Nj) = (2048, 2048, 32), 
(1024,1024, 32), (1024,1024, 512) and (512, 512, 512) from top to bottom, as indicated on each panel. The pictures show only the / ^ 0 
part of the phase-space density, while it can actually become negative because of aliasing. However, this choice of representation does 
not hide aliased regions. The prominent one corresponds to the textured zone above the large r tail of the system on the right panels. 
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Figure A3. Same as in Fig. |A2l but the phase-space distribution function has now been summed up over the whole available range of 
values of j G [0, Jmax = 1-6], where Jmax is the maximum sampled value of j. 
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Figure Bl. Effect of changing the important control parameters in Gadget. The phase-space density is shown at t = 15 for Gadget 
simulations with the same initial conditions corresponding to the Henon sphere with = 0.1 and involving N = 10® particles. In each 
of the simulations, one control parameter was changed compared to the fiduciary simulation shown on left panel and which uses the 
settings of § [^ On top and bottom left panels, the softening length of the force was decreased by a factor 5 and increased by a factor 
10, respectively. In top-right panel, the maximum possible time step was divided by a factor 50, while in the bottom-right panel, the 
tolerance parameter ap defined in §[^was divided by a factor 5. 
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APPENDIX B: A-BODY SIMULATIONS: 
EXPLORATION OF THE CONTROL 
PARAMETER SPACE 

In §[^we noticed the presence of an instability in the R = 0.1 
iV-body simulations. One aim of this appendix is to confirm 
that this instability is related to the number of particles used 
in the simulations and not to any other control parameter 
of the Gadget code. In the same time, it is also an opportu¬ 
nity to check that our fiducial choice of the Gadget control 
parameters, given in §[^ is correct. 

Figure |B1| illustrates the main results of the tests we 
performed for simulations with 10® particles. These tests 
consisted in changing the softening length of the force, the 
maximum time step value and the tolerance parameter ap 
controlling the errors on the force. Improving the accu¬ 
racy of the force calculation or dividing the maximum time 
step dtmax by a factor 50, which corresponds to imposing 
dt ^ 2 X 10“^, does not change the results. This is con¬ 
firmed as well by the measurements of the correlators Ck 
introduced in § [^ that we do not show here for simplicity. 
Only the value of the softening parameter of the force e has 
an impact on the dynamics for the tests we did. Reducing e 
by a factor 5 seems to slightly blur the phase-space density, 
although this effect is difficult to decipher, while increasing e 
by a factor 10 sharpens the fine structures of the phase-space 
density. Since e controls the intensity of close encounters be¬ 
tween particles, this is not surprising. Note that increasing e 
by a factor 10 is probably an exaggeration, because it wors¬ 
ens dramatically the match during the mixing phase between 
the A-body simulation and the Vlasov code when examin¬ 
ing the correlators C/c, a sign that e is probably getting too 
close to a physical characteristic scale of the systemjjWe 
indeed noticed that increasing e only by a factor 5 does not 
have much impact, on the other hand, on Ck- However, all 
these effects do not affect the amplitude of the large scale 
irregularities on the pattern of f{r,v,j), which are present 
whatever value of e. This is also a strong indication that 
close particle encounters are not at the origin of these irreg¬ 
ularities. 

We can therefore only conclude that these irregularities 
and the associated nonlinear instability are the result of non 
trivial collective effects related to particle shot noise. This 
argument is also supported by the fact that in addition, the 
moment of their appearance is particle number dependent, 
as discussed in § [^ 


^ Increasing e by a factor ten gives e = 0.02, to be compared 
for example to the size of the core of the system after relaxation, 
Rc ^0.1. 
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